Core sizes and dynamical instabilities of giant vortices in dilute Bose-Einstein 

condensates 
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Motivated by a recent demonstration of cyclic addition of quantized vorticity into a Bose-Einstein 
condensate, the vortex pump, we study dynamical instabilities and core sizes of giant vortices. The 
core size is found to increase roughly as a square-root function of the quantum number of the vortex, 
whereas the strength of the dynamical instability either saturates to a fairly low value or increases 
extremely slowly for large quantum numbers. Our studies suggest that giant vortices of very high 
angular momenta may be achieved by gradually increasing the operation frequency of the vortex 
pump. 

PACS numbers: 03.75.Kk, 03.75.Lm, 67.85.De 



I. INTRODUCTION 

One of the signatures of Bose-Einstein condensation is 
the occurrence of superfluidity. In superfluids, particle 
currents can flow without dissipation, and the fluid is 
characterized by zero viscosity. A stable quantized vor- 
tex is a hallmark of such frictionless flow. These coherent 
whirlpools manifest the fundamental differences in the 
rotational characteristics of quantum and classical fluids. 
The study of quantized vortices [l[ has flourished espe- 
cially after their first realization in dilute Bose-Einstein 
condensates (BECs) Q. 

Large arrays of singly quantized vortices have been 
observed in rapidly rotating clouds of cold alkali metal 
atoms Furthermore, multiquantum vortices, for 

which the phase of the condensate order parameter winds 
an integer multiple k of 2tt, have been created in rotat- 
ing clouds by using an evaporative cooling technique to 
spin up the condensate Q in combination with a resonant 
laser beam focused at the center of the cloud § . The gi- 
ant vortices were reported to have surprisingly long life- 
times owing to strong Coriolis forces. However, the con- 
densate states obtained in these experiments were not 
close to a pure, symmetric multiquantum vortex state, 
since the giant vortices were accompanied by several 
single-quantum vortices. Dynamics of such giant vortex 
states have since been investigated both experimentally 
and theoretically d, [l(| . 

Multiquantum vortices have also been created in an 



initially nonrotating cloud by transferring angular mo- 
mentum into the condensate by a Laguerre-Gaussian 
laser bea m fli"l .[l2| and by a topological phase engineering 
method [13[ which utilizes the spin degree of freedom of 
the condensate and its coupling to an external magnetic 
field [31(1. States with a multiquantum vortex are in 
general dynamically unstable: the state may annihilate 
due to a slight perturbation even in the absence of dis- 
sipation. Dynamical stability of multiquantum vortices 
has been investigated theoretically j2l!425} . and splitting 
of multiquantum vortices into single-quantum ones has 
been studied both experimentally and theoretically [26l — 
. In addition, dynamical stability of coreless vor- 
tices [H, HI1 and vortex clusters (35l - l40j has been inves- 
tigated. 

Lately, it has been suggested that vortices with arbi- 
trarily large winding numbers could be created by cycli- 
cally pumping vorticity into a BEC [4l[ . The technique is 
based on the topological phase engineering method: the 
spin degree of freedom of the condensate is controlled 
locally by alternating external quadrupole and hexapole 
magnetic fields, and a fixed amount of vorticity is added 
into the system in each cycle. In the simulations of 
Ref. " 
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the operation of the vortex pump was demon- 
strated both adiabatically, in which case an additional 
optical potential is required to confine the condensate 
during the pumping cycle, and nonadiabatically with- 
out any optical potentials, but with the cost of losing 
a part of the atoms from the trap. Later, the vortex 
pump was investigated numerically using only a single 
multipole magnetic field, albeit an additional uniform 
field was introduced [42} ■ This scheme has the advantage 
that no added vorticity is removed during the cycle, but 
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it necessarily requires an optical potential alongside the 
magnetic fields. Recently, it has also been theoretically 
demonstrated that topological phase engineering can be 
used to create a Dirac monopole in a spinor BEC (43j . 

In this paper, we study the vortex core sizes and dy- 
namical instabilities of axisymmetric giant vortex states 
in dilute nonrotated BECs as functions of the vorticity 
quantum number and the effective atom-atom interac- 
tion strength. We restrict our investigation to the zero- 
temperature limit and neglect the possible additional ef- 
fects due to thermal atoms. Hence, the analysis presented 
here provides an upper bound for the lifetime of the giant 
vortices. 

Our investigation is motivated by two aspects. Firstly, 
the studies of multiquantum vortices in dilute BECs to 
date have focused on the regime K < 5, probably be- 
cause previously there has been no practically realizable 
techniques available to create isolated vortices with large 
quantum numbers in a controlled manner. Now, how- 
ever, the vortex pump provides such an opportunity. Sec- 
ondly, we aim to determine how large winding numbers 
the pump can reach for a single vortex. In principle, 
vortex pumping enables a controlled production of vor- 
tices with arbitrarily large quantum numbers. However, 
in real experiments the quantum number will be limited 
by the stability properties of the vortex states. In par- 
ticular, the existence of dynamical instabilities renders 
the multiquantum vortices prone to splitting into single- 
quantum vortices. The frequencies of these instabilities 
can be used to approximate the lifetimes of the multi- 
quantum vortex states. On the other hand, the core size 
of the vortex is known to increase with increasing wind- 
ing number. As discussed in Sec. IIV1 this enables us 
to speed up the pumping cycle as the vorticity of the 
condensate increases. Hence, the interplay between the 
core size and the lifetime of the vortex as a function of 
the winding number ultimately determines the maximum 
winding number reachable with the vortex pump. 

The remainder of this paper is organized as follows. 
We present our theoretical model in Sec. Inland also con- 
sider the Thomas-Fermi (TF) approximation for giant 
vortex states. In Sec. IIII1 we formulate the analytical 
approximation in the noninteracting limit. Section [IVI is 
devoted to presenting our numerical results on the core 
sizes and dynamical instabilities of multiquantum vor- 
tices. We compare the numerically obtained core sizes 
with the analytical results of the TF approximation and 
consider how well the noninteracting limit can be used 
to describe the dynamical instabilities in the presence of 
finite atom-atom interactions. In Sec. [Vj we summarize 
and discuss the main results of our study. 



II. THEORETICAL FRAMEWORK 

In the zero-temperature limit, the complex-valued or- 
der parameter field ^(r, t) of a dilute BEC is described 



by the time-dependent Gross-Pitacvskii (CP) equation, 



5 |*M)| 2 ]vE-(r,i), 



(1) 



where % denotes the single-particle Hamiltonian and the 
atom-atom interaction strength parameter g is related to 
the vacuum s-wave scattering length a and the atomic 
mass to by g — Airh 2 a/m [4J|. For BECs enclosed in 
nonrotated cylindrically symmetric harmonic traps, the 
Hamiltonian % is expressed in the cylindrical coordinates 
(r, </>, z) as 



2m 



V 2 + -tow 2 (r 2 + AV) , 



(2) 



where uj r is the harmonic oscillator frequency in the ra- 
dial direction and A = uj z /uj r defines the aspect ratio of 
the trap. The order parameter is normalized according 
to J \^(r)\ 2 dr = TV, with N being the total number of 
particles. 

Stationary states of the system are solutions to the 
time-independent CP equation, 



[H + .g|v]/(r)| 2 - M ]*(r) = 0, 



(3) 



where fi denotes the chemical potential. In order to 
study the small-amplitude oscillations about a given sta- 
tionary state, we write the order parameter in the form 
^(r, t) = exp(—ifj,t/h) [^(r) + %(r, t)} , where x(r, t) is as- 
sumed to have a small L 2 -norm compared with ^(r). By 
substituting this trial function into Eq. ([Tj), using the 
decomposition 



X(r.t) 
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\u q {r)t 



v*(r)e 
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(4) 



and linearizing with respect to x(r,i), we arrive at the 
Bogoliubov equations 
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u q {r) 
v q (r) 



(5) 



where we have denoted C = H + 2g\^f\ 2 — /i. An al- 
ternative derivation utilizing a canonical transformation 
of second-quantized operators shows that the functions 
u q (r) and v q (r) are the quasiparticle amplitudes and u> q 
is the eigenfrequency of the elementary excitation mode 
corresponding to the index q (44j . 

The Bogoliubov equations can be used for the sta- 
bility analysis of the system. If the quasiparticle spec- 
trum {uj q } contains excitations with a positive norm 



J[N S 



2 ] dr but a negative eigenfrequency ui q , the 



corresponding stationary state is energetically unstable. 
On the other hand, the stationary state is dynamically 
unstable if the quasiparticle spectrum contains at least 
one eigenfrequency with a positive imaginary part [2l| . 
As can be observed from Eq. (|4]) , the occupations of such 
complex-frequency modes initially increase exponentially 
in time, and consequently small perturbations of a dy- 
namically unstable stationary state typically lead to large 
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changes in its structure. In particular, dynamically un- 
stable multiquantum vortex states are unstable against 
splitting of the vortex, and for such states the quantity 
max g [|Im(w 9 )|] can be used to estimate the splitting ten- 
dency of the vortex [27]] • However, the dynamically un- 
stable modes quickly drive the system beyond the linear 
regime of Eq. Q , and hence the dynamics must instead 
be described with the time-dependent GP equation. 

In this paper, we study the properties of stationary 
multiquantum vortex states as functions of the winding 
number k of the vortex. Since the external potential 
in Eq. ([5]) is rotationally symmetric, a stationary state 
containing a K-quantum vortex can be written in the form 



(6) 



where n(r, z) is the particle density of the condensate. 
In order to simplify the analysis, we limit our consid- 
erations to pancake-shaped BECs, i.e., we assume that 
the harmonic confinement in the axial direction is suf- 
ficiently tight such that the axial harmonic oscillator 
length a z = hjmuiz is much smaller than the char- 
acteristic length scale of density variations in the radial 
direction. This approximation is accurate for A> 1, and 
in such a case the z dependence of the order parame- 
ter can be factored out as y/ n(r, z) = f(r)(o(z), where 
Co(^) = cxp[— z 2 / (2a 2 z )]/ y/ira 2 . Moreover, for this kind 
of stationary states, the azimuthal and axial dependences 
can be separated from the solutions of the Bogoliubov 
equations ([5]), and consequently the quasiparticlc ampli- 
tudes can be expressed in the form 

u g (r,4>,z) = u q (r)Co(z)e^ + ^, 

v q (r,</>,z) = v q (r)C (z)e^- K ^, (7) 

where l q is an integer that determines the angular mo- 
mentum of the quasiparticle excitation with respect to 
the condensate. 

It is convenient to express all quantities using dimen- 
sionless variables, which are henceforth denoted with a 
tilde. Therefore, we measure leng th in th e units of the 
harmonic oscillator length a r = y/h/mujr and energy in 
units of hiu r . Using the factorization of Eq. © with 
yj n(r, z) = f(r)Ca(z), we arrive at the one-dimensional 
GP equation 



n K + ~gf\f)-fi /(f) = 0, 



(8) 



bility density of the condensed particles. The dimension- 
less Bogoliubov equations for a ^-quantum vortex state 
are given by 



£i 1+ k gf 2 (?) 
-gf 2 (r) -Cu- K 



Uq(f) 



Uq(f) 
Vqif) 



(10) 



where C K = H K + 2gf 2 (r) — jl. Using the values 
a = 2.75 nm and m = 3.81 x 10 -26 kg corresponding 
to 23 Na atoms [44| and taking the radial oscillator fre- 
quency and the trap asymmetry to be uj r — 2ix x 20 Hz 
and A = 50, which are typical values for traps used in ex- 
periments [45[ , we obtain a r « 4.7 /^m and g tn 0.021 x N. 
Thus, for a 23 Na condensate with N ps 10 4 atoms, the 
effective interaction strength is of the order of g 200. 

For sufficiently large numbers of atoms, i.e., for large 
enough values of <?, a plausible description for the con- 
densate order parameter is obtained by solving the GP 
equation in the Thomas-Fermi (TF) approximation [46| . 
In our case, the TF solution corresponds to neglecting 
the radial derivatives of /(f) in Eq. ((5J), which yields 



n TF (f) 



K 



i 



Mtf 



f 2 



lK 2 

2f 2 



(11) 



if ftTF(f) > 0, and otherwise. The inner and outer TF 
radii, between which nTp(f) is nonvanishing, are given 
by f% = Atf T \J Mtf — k2 - The TF expression for the 
chemical potential, /UtFi is determined by the normal- 
ization condition 2-7T J nTF(f)fdf = 1, which yields the 
equation 



.9 _ ~2 
— — Mtf 

7T 



1 - 



hi 



Mtf^ 4 , (12) 



where 9 = y/l - k 2 /^. In Eq. (JUJ, we have approx- 
imated the bracketed function by 9 , which is accurate 
for l/\/2 < 9 < 1 (the relative error is less than 7%). 
Finally, solving for the chemical potential gives 



9 I 9 

MTF = \/^ + J^ 



(13) 



and thus the approximation in Eq. (1121) is justified if g > 

k 2 m. 



where the dimensionless interaction strength parame- 
ter is defined as g = 2y/2~nNa/a z , the chemical poten- 
tial jl = ^,/{huj r ) — ui z /(2uj r ) contains the shift due to 
the z displacement, and the dimensionless single-particle 
Hamiltonian is given by 



1 fld_ 

2 \ f dr 



dr 2 



f 2 



(9) 



III. NONINTERACTING LIMIT 

In the limit of weak interactions, the instabilities can 
be calculated by treatin g th e interaction strength g as a 
perturbative parameter |25l . |48| . The unperturbed wave- 
functions are the single-particle eigenfunctions without 
radial nodes, 



Function / is normalized according to the condition 
2ir J h(r)rdf = 1, where n(f) = / 2 (f) is the areal proba- 



<Pn(r) 



:rl"le- f / 2 , 



(14) 
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where n £ Z and the angular dependence exp(m</>) has 
been factored out. Thus, for a K-quantum vortex state, 
we haw 1 



(15) 



and the corresponding Bogoliubov eigenstates are the 
particle states 



Uq(f) \ _ ( tp K +i q (r) 

«5(r) 



(16) 



with energies tb q = \k + l q \ — k, and the hole states 



vPJf) 
v\{f) 



(17) 



with energies ui q = — 1« — l g \ + K. We note that for l q < k, 
the particle and hole states with equal angular momen- 
tum quantum numbers l q are degenerate. By applying 
degenerate-state perturbation theory, we find the Bogoli- 
ubov energies 



IV. RESULTS 

We have solved the stationary states and the cor- 
responding Bogoliubov excitation spectra, Eqs. © 
and (|10p. for the parameter ranges < k < 100 and 
< g < 10 4 . In the calculations, we have used finite 
difference methods with grid sizes of roughly 200 points 
and LAPACK numerical library for solving the discretized 
eigenvalue problem. 

In order to illustrate a typical profile of a giant vortex 
state and its TF approximation, Fig. [TJ shows the mod- 
ulus of the order parameter for g = 1000 and different 
values of k, together with the corresponding TF profiles. 
As expected, the TF approximation fails near the inner 
and outer surfaces of the condensate. Inside the cloud, 
the TF profile provides an accurate estimate as long as 
the condition g > k 2 is satisfied. The state with n = 40 
does not satisfy this criterion, and consequently the TF 
profile significantly deviates from the numerical solution 
also away from the surface of the cloud. 



Size of the vortex core 



"9 = k + .9^7-^- ± W ( J + + - ( 18 ) 



where 



I± = 2ir J <f K (?) 2 cp K±lq (?) 2 rdr, (19) 
h = 2tt J ^ K (f) 4 fdf, (20) 



and 



I x = 2tt / (p K {r) v K+lq (r)<p K -i g (r)rdr. (21) 

The integrals can be evaluated using the formula [48| 

2n f <p m (r)<p n (r)<Pk(r)ipi(r)r dr = 
1 [(|m| + |n| + |fc| + |l|)/2]! 



2tt 2(l m l+l"l+l fe l+l'IV2 A /| TO |!| n |!|A;|!|z|!' 



(22) 



We are interested in the maximum imaginary part of 
the Bogoliubov eigenfrequencies. In the noninteracting 
limit, it is found by minimizing the expression inside the 
square root in Eq. (|T5|) . Note that the perturbation the- 
ory does not describe the eventual disappearance of a 
given unstable mode at higher values of g, but incor- 
rectly predicts that the instabilities grow linearly with 
g. Nevertheless, for a given winding number k, one can 
determine the angular momentum quantum number l\ om 
that gives the dominant instability in the noninteracting 
limit. These analytical results are presented and com- 
pared to the numerical ones in Sec. IIVI 



Here, we study the radius of the vortex core, r c , as 
a function of the winding number k and the interaction 
strength g. We define the radius as the smallest solution 
of the equation 



f 2 (r c ) = -max/ 2 (f) 

4 r 



(23) 



where f 2 (f ) is the areal probability density of the station- 
ary K-quantum vortex state satisfying Eq. (jHJ . Our inves- 
tigation of the vortex core size is partly motivated by the 
fact that for the vortex pump [4l| , the core radius essen- 
tially determines the maximum adiabatic pumping speed 
in the following way. A characteristic energy associated 
with the magnetic field is given by the minimum energy 
separation e mag = gFfJ-BB m i n , where g-p is the Landc g 
factor, //b is the Bohr magneton, and B nlln is the min- 
imum magnetic field strength in the condensate region 
during the bias field inversion. The field strength B is 
minimized when the spatially homogeneous bias field B z 
crosses zero, at which point the magnitude of the field 
increases linearly with distance from the z axis. Thus, 
the relevant minimum field strength in the presence of 



vortex becomes B„ 



B 



_L'c, 



where B' ± denotes the 



gradient of the perpendicular multipole field. To guaran- 
tee adiabaticity, the time derivative of the Hamiltonian 
should be small compared with the energy separation of 
its instantaneous eigenstates. In our case, this roughly 
means that the bias field 6, should be inverted such that 



ftffFMB dB z 
— n x 



'-mag 



at 



g F fj, B (B' ± ) 



3B Z 



(24) 



Since the core radius r c is observed to increase with k, it 
is possible to speed up the pumping cycle without chal- 
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longing adiabaticity as more vortices are pumped into the 
BEC. 

In order to obtain an analytical approximation for the 
vortex core size, we use the TF density profile, Eq. (1111) . 
Inserting Eq. (fTTj) into Eq. (j2"5)) . we find the expression 



rJ F = - y Atf + 3k - y /t| F + 6k/2tf - 7k 2 , (25) 

which, with the help of Eq. (fT3"|). can be directly com- 
pared with the numerical results. Furthermore, since the 
TF inner radius f_ gives an estimate for the radius be- 
low which there are no particles, it could also be used to 
evaluate the adiabaticity condition of the vortex pump 
as discussed above. 

On the other hand, the core size of single-quantum 
vortices is typically approximated by the healing length 
£ = l/\/87ma, which describes the characteristic distance 
over which the condensate density tends to its bulk value 
n when subjected to a localized perturbation in an oth- 
erwise homogeneous system [44j . For multiquantum vor- 
tices, we can use a similar approach to obtain a simple 
approximation for the core radius. However, since the 
size of the giant vortex is comparable to that of the whole 
condensate, the inhomogeneity of the trapped gas can no 
longer be neglected, but instead the trapping potential 
term must be included in the energy balance. Hence, we 
equate the kinetic energy with the potential and interac- 
tion energy terms at f = £ c and obtain 

■i = + U, (26) 
2^ c v ' 

where U denotes the interparticle interaction. For a 
vortex-free condensate in a harmonic potential in the 
TF limit, one finds that near the center of the cloud 
U oc TV 2 / 5 (44[. Since g oc N, we make the crude es- 
timate U ~ g 2 / 5 to obtain the core radius approximation 

| c = vV« 2 +.9 4 / 5 -.9^, (27) 

which should be compared with the more rigorous ap- 
proximation given by Eqs. (|2"5"j) and (fT3")) . We note that 
for large quantum numbers k ^> g 2 ^ 5 , Eqs. (f2"5j) and (|27|) 
predict that the core radius increases as a square root of 

K. 

In Fig. [2] we plot the vortex core radius as a function 
of k for g — 1 and 10 4 together with the approximations 
of Eqs. ([231) an d (|27|) . Figure |3] shows the radius versus 
the interaction strength g for different winding numbers. 
As expected, f c monotonously increases with increasing 
k and decreases with increasing g. As a function of n, 
the asymptotic behavior f c cx *Jk is observed for k 2 ^> g, 
as predicted by Eqs. (f2"5j) and (|2T| . The reduction of the 
core size with the interaction strength is relatively slow, 
and f c decreases by a few tens of percent over the range 
< g < 10 4 . In general, the analytical TF expression for 
the vortex core radius, Eq. (f2"5"]). yields good agreement 
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FIG. 1: (Color online) Numerical solution (solid curve) and 
the TF approximation (dashed curve) for /(f) as a function 
of f for g — 1000. The curves correspond to the winding 
numbers n = 1, 10, 20, 40 as indicated. 
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FIG. 2: (Color online) Radius of the vortex core as a func- 
tion of the winding number k for g = 1 and g — 10 4 . The 
solid curve shows the numerically obtained radius f c , the dot- 
ted curve is the TF approximation f J F given by Eqs. (|13|) 
and (|25[) , and the dashed curve corresponds to the approxi- 
mation £ c , Eq. (|27[) . For g = 10 4 , the solid curve lies on top 
of the dotted curve. 

with the numerical results, even for small values of g. On 
the other hand, the simpler approximation of Eq. (|27[) 
is observed to give accurate results for sufficiently small 
values of g and large values of k. 



B. Dynamical instabilities 

Next, we investigate the dynamical instability proper- 
ties of the vortex states as functions of the winding num- 
ber and the interaction strength. It is known from ear- 
lier numerical simulations that in pancake-shaped con- 
densates, the multiquantum vortex can be dynamically 
stable or unstable depending on the interaction strength 
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FIG. 3: (Color online) Radius of the vortex core as a func- 
tion of the interaction strength g. The solid curve shows 
the numerically obtained radius f c , the dotted curve is the 
TF approximation fJ F , and the dashed curve corresponds to 
the approximation £ c . The curves correspond to the winding 
numbers k = 10, 30, 50, 100 as indicated. 



9 [2ll. |23||. In a harmonic trap, the two-quantum vortex 
state can have a dynamical instability mode only for the 
angular momentum quantum number \l q \ = 2. A vor- 
tex with a higher winding number k can have complex- 
frequency modes also at larger values of \l q \, each mode 
corresponding to an |Z g |-fold symmetric splitting pat- 
tern 
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These observations are also predicted by the 
noninteracting approximation of Sec. IIII1 according to 
which degenerate modes exist for those l q that satisfy 
l q < k. Because the number of different Z 9 -modes sup- 
porting complex cigenfrequencies tends to increase with 
K, we expect the multiquantum vortices to generally be- 
come more and more unstable as the vorticity increases. 

Since we are interested in the general degree of dy- 
namical instability of a vortex state with a given wind- 
ing number k, we need a quantity that is independent 
of the interaction strength g. Therefore, we define the 
dimensionlcss instability parameter T K by 



T K = max < max [|Im( 



(28) 



In other words, we solve Eq. (fTU)) for fixed values of k and 
g and maximize the quantity |Im(o) g )| over each eigen- 
frequency spectrum, after which we vary the interaction 
strength and maximize max,j[|Im(w g )|] as a function of g. 
The parameter T K provides a measure for the splitting 
tendency of a K-quantum vortex, and 1 /T K can be used 
to estimate the lifetime of the vortex state. Furthermore, 
it was found in Ref. [13] that a highly elongated BEC, 
i.e., one with a small aspect ratio A, can be described by 
a local-density approximation in which the condensate is 
viewed as a stack of two-dimensional systems with differ- 
ent values of g. If the conditions for two-dimensional in- 
stability are met locally at a point z, the vortex splitting 
will commence at that point and subsequently migrate 
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FIG. 4: (Color online) Dimensionless dynamical instability 
parameter T K of a giant vortex state as a function of its 
winding number k. The dashed vertical lines separate re- 
gions where the angular momentum quantum numbers of the 
dominant complex- frequency modes, idom, are different. 
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FIG. 5: (Color online) Values of the dimensionless interaction 
strength g at which the largest dynamical instabilities occur. 
The angular momentum quantum number of the dominant 
complex- frequency mode, Idom, is also indicated. 



along the vortex line. Thus, finding the maximum of the 
quantity maxg[|lm(a; 9 )|] with respect to g is equivalent 
to finding the dominant unstable mode of a cigar-shaped 
BEC in this approximation. 

The numerically obtained values of T K are presented 
in Fig. |U Figure [5] shows the values of the interac- 
tion strength, <7dom, for which the dominant complex- 
frequency mode occurs. The angular momentum quan- 
tum number /dom of the dominant mode is also indicated. 
As expected, T K is a strictly increasing function of k. 
However, we observe that the value of T K cither saturates 
to a relatively low value of about 0.3 or increases very 
slowly for large winding numbers. In addition, the in- 
stability parameter exhibits quasiperiodic behavior that 
results from the stepwise increments in the value of Zdom- 
Quasiperiodic behavior is observed also for c/dom in Fig. [S] 
<?dom increases within regions of constant idom but ex- 
hibits a sudden drop each time Idom increases. 
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FIG. 6: (Color online) Dots show the angular momentum 
quantum numbers /^ om yielding the largest dynamical insta- 
bility in the noninteracting limit. The triangles denote the 
critical interaction strength, g°r, above which the noninter- 
acting result for idom fails. The dashed line is a guide to the 
eye. 



In earlier studies for 2 < k < 5, the dominant complex- 
frequency mode over different values of the interaction 
strength was found at the angular momentum quantum 
number l q = 2. According to our results, the l q = 2 mode 
dominates up to the winding number k = 20, after which 
the largest imaginary part is provided by an excitation 
mode with l q = 3. The subsequent transition to /dom = 4 
occurs at k = 70. These findings should be compared 
with the analytical results, displayed in Fig. [51 for the 
angular momentum quantum number lJ om that provides 
the largest complex-frequency mode in the noninteract- 
ing limit for a given winding number. We have also indi- 
cated the numerically obtained value of g above which the 
angular momentum quantum number of the most unsta- 
ble mode becomes different from l% om - We observe that 
!j om increases slowly as a function of k. For small k, 
the most unstable mode has the quantum number l^ om 
for exactly l° om consecutive values of K. However, this 
numerical coincidence ceases to hold after l° om = 7. We 
observe that in the noninteracting limit, the dominant 
complex-frequency modes generally have much larger an- 
gular momentum quantum numbers than is found numer- 
ically for unrestricted g. Moreover, the range of validity 
for the noninteracting result is restricted to very small 
values of g. 



V. DISCUSSION 

In conclusion, we have studied the core sizes and 
dynamical instabilities of multiquantum vortices in 



pancake-shaped nonrotatcd BECs as a function of the 
winding number n. We found the core size to be a 
strictly increasing function of the winding number with 
the asymptotic behavior f c ex ^/k for large values of 
k. The dependence of r c on the interaction strength g 
turned out to be weak. On the other hand, the maximum 
strength of the dynamical instability of a K-quantum vor- 
tex state was observed to increase slowly with k or even 
to saturate to a relatively low value. Based on these re- 
sults, the time scales for the splitting of giant vortices are 
expected to decrease only slowly with increasing winding 
number. In addition, we found that quasiparticle exci- 
tations with the angular momentum quantum number 



L 



2 dominate the instability spectra up to the wind- 



ing number n — 20, after which consequent transitions 
to higher modes are observed. 

We found that the analytical TF approximation pro- 
vides good estimates for the exact particle densities as 
long as the condition g > k 2 is satisfied. Moreover, the 
TF estimate for the radius of the multiquantum vortex 
was found to be in good agreement with the numerically 
determined radius for all parameter values {k, g} consid- 
ered. 

From the point of view of the vortex pump (4l| , our re- 
sults are encouraging. Qualitatively, the efficiency of the 
pump in producing vortices with large winding numbers 
is described by the ratio T sp /Ti nv , where r sp is a char- 
acteristic time scale for vortex splitting and T[ nv is the 
bias field inversion time. Since the radius of the vortex 
core increases with the winding number approximately as 
r c ex it should be possible to gradually decrease T; nv 
and still retain the adiabaticity of the process. On the 
other hand, we expect the splitting times r sp to decrease 
only slowly with k for large vorticitics. Consequently, it 
should be possible to actually increase the ratio r sp /Ti nv 
after a sufficient amount of vorticity has accumulated 
into the condensate. Hence, our results suggest that gi- 
ant vortices with very high quantum numbers may be 
achieved by gradually increasing the operation frequency 
of the pump, e.g., such that T; nv oc 
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